library("rstudioapi")     
setwd(dirname(getActiveDocumentContext()$path))
#setwd("../")
getwd()
library("readxl") 
require(miceadds)
library(data.table)
library(tidyverse)
library(DescTools)
library(htmlTable)
library(stargazer)
library(xlsx)
library(estimatr)
library(tidyverse)
library(DescTools)
require(dplyr)

load("FinalFinal20062023Full1-6.RData")
load("FinalFinalSP20062023Full1-6.RData")

library(dplyr)
final_finalV2 = mutate(final_finalV2, 
                       Agec = (Age - 30)/10, # per 10 year age
                       Incomec = (Income - 100)/50, # per 50 income
                       clinic_distancec = (clinic_distance - 5)/5, # per 5
                       SNMetricc = SNMetric/10)

L_7c <- glm.cluster(vaccine_intention ~ cash_cat + Gender + Agec + 
                      Education + Incomec + Employed + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) 

L_8c <- glm.cluster(vaccine_reported_combo ~ cash_cat  + Gender + Agec + 
                      Education + Incomec + Employed + SNMetricc + clinic_distancec +
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) 

L_9c <- glm.cluster(ActVacApril ~ cash_cat  + Gender + Agec + 
                      Education + Incomec + Employed + SNMetricc + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123)


tab_7c1 = exp(cbind("Odds ratio" = coef(L_7c), confint.default(L_7c, level = 0.95)))
tab_8c1 = exp(cbind("Odds ratio" = coef(L_8c), confint.default(L_8c, level = 0.95)))
tab_9c1 = exp(cbind("Odds ratio" = coef(L_9c), confint.default(L_9c, level = 0.95)))
rownames(tab_7c1) = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Access","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_8c1) = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_9c1) = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access","District 2", "District 3", "District 4", "District 5","District 6")

# Convert function, to create outputs for stargazer tables: 
convert_fun = function(data, name_mod){
  new_vec = c()
  for (i in c(1:dim(data)[1])){
    new_vec[i] = as.character(paste(format(round(data[i,1],2), nsmall = 2),paste("(",paste(format(round(data[i,2],2), nsmall = 2), format(round(data[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
  }
  new_vec = data.frame(new_vec)
  rownames(new_vec) = rownames(data)
  colnames(new_vec) = name_mod
  return(new_vec)
}

col_7c1 = convert_fun(tab_7c1, "Model1")
col_8c1 = convert_fun(tab_8c1, "Model1")
col_9c1 = convert_fun(tab_9c1, "Model1")

# Other models: 
### Ghana - Attrition: 
final_finalV2$full=ifelse(final_finalV2$Q143 == "Employed (full time)", 1 , 0)
final_finalV2$part_time=ifelse(final_finalV2$Q143 == "Employed (part time)", 1 , 0)
final_finalV2$unemploy=ifelse(final_finalV2$Q143 == "Unemployed", 1 , 0)

final_finalV2$p_placebo=ifelse(final_finalV2$individual_treatment == "Placebo", 1, 0)
final_finalV2$p_health=ifelse(final_finalV2$individual_treatment == "Placebo", 0 ,
                              ifelse(final_finalV2$individual_treatment == "CDC Health", 1 , 0))
final_finalV2$p_low=ifelse(final_finalV2$individual_treatment == "Placebo", 0 ,
                           ifelse(final_finalV2$individual_treatment == "Low Cash", 1 , 0))
final_finalV2$p_high=ifelse(final_finalV2$individual_treatment == "Placebo", 0 ,
                            ifelse(final_finalV2$individual_treatment == "High Cash", 1 , 0))

final_finalV2$id <- 1:nrow(final_finalV2)
final_finalV2$cens <- as.numeric(is.na(final_finalV2$vaccine_reported_combo)) # turns TRUEs into 1 and 0 if not missing
table(final_finalV2$cens)

library(nnet)
### Model combining all treatments (categories 1/0 for each treatment): 

denom.cens <- glm(cens ~ 
                    p_health:Age + p_health:Gender + p_health:Education + 
                    p_low:Age + p_low:Gender + p_low:Education +
                    p_high:Age + p_high:Gender + p_high:Education,
                  data=final_finalV2, family=binomial)
summary(denom.cens)

denom.p.cens <- predict(denom.cens, type = "response") # gets the predicted probabilities 
final_finalV2$wt <- ifelse(final_finalV2$cens == 0, (1/(1 - denom.p.cens)),1) 
# for not missing values, if gets the predicted probabilities, for missing values it gives 1

summary(final_finalV2$wt[final_finalV2$p_placebo=="0"])
summary(final_finalV2$wt[final_finalV2$p_placebo=="1"])

# re-estimate models: 
final_finalV2$sel_valid <- final_finalV2$cens == 0 

Mod1 <- glm(vaccine_intention ~ cash_cat + Gender + Agec + 
              Education + Incomec + Employed + clinic_distancec+
              dist2 + dist3 + dist4 + dist5 + dist6, data = final_finalV2, 
            weights = wt, subset=sel_valid)
summary(Mod1)

Mod2 <- glm(vaccine_reported_combo ~ cash_cat + Gender + Agec + 
              Education + Incomec + Employed + SNMetricc+ clinic_distancec+
              dist2 + dist3 + dist4 + dist5 + dist6, data = final_finalV2,
            weights = wt, subset=sel_valid)
summary(Mod2)

Mod3 <- glm(ActVacApril ~ cash_cat + Gender + Agec + 
              Education + Incomec + Employed + SNMetricc+ clinic_distancec+
              dist2 + dist3 + dist4 + dist5 + dist6, 
            data = final_finalV2,
            weights = wt, subset=sel_valid)
summary(Mod3)

mm1 = exp(cbind("Odds ratio" = coef(Mod1), confint.default(Mod1, level = 0.95)))
mm2 = exp(cbind("Odds ratio" = coef(Mod2), confint.default(Mod2, level = 0.95)))
mm3 = exp(cbind("Odds ratio" = coef(Mod3), confint.default(Mod3, level = 0.95))) 

rownames(mm1) = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Access","District 2", "District 3", "District 4", "District 5","District 6")
rownames(mm2) = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access","District 2", "District 3", "District 4", "District 5","District 6")
rownames(mm3) = c("Intercept", "Cash", "Health", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access","District 2", "District 3", "District 4", "District 5","District 6")


col_m1 = convert_fun(mm1, "Model1")
col_m2 = convert_fun(mm2, "Model1")
col_m3 = convert_fun(mm3, "Model1")


tb1 = merge(col_7c1, col_m1, by = 0, all = TRUE);rownames(tb1) = tb1[,1];tb1[,1] = NULL
tb2 = merge(col_8c1, col_m2, by = 0, all = TRUE);rownames(tb2) = tb2[,1];tb2[,1] = NULL
tb3 = merge(col_9c1, col_m3, by = 0, all = TRUE);rownames(tb3) = tb3[,1];tb3[,1] = NULL
ttb1 = merge(tb1, tb2, by = 0, all = TRUE);rownames(ttb1) = ttb1[,1];ttb1[,1] = NULL
ttb2 = merge(ttb1, tb3, by = 0, all = TRUE);rownames(ttb2) = ttb2[,1];ttb2[,1] = NULL
colnames(ttb2) = c("7c", "7c_weight", "8c", "8c_weigth", "9c", "9c_weight")
no.obs = c((L_7c$glm_res$df.null),(Mod1$df.null), (L_8c$glm_res$df.null), (Mod2$df.null), (L_9c$glm_res$df.null), (Mod3$df.null))
aic_mod = round(c(L_7c$glm_res$aic, Mod1$aic, L_8c$glm_res$aic, Mod2$aic, L_9c$glm_res$aic, Mod3$aic),2)
log_lik = round(c(as.vector(logLik(L_7c$glm_res)), as.vector(logLik(Mod1)), as.vector(logLik(L_8c$glm_res)), as.vector(logLik(Mod2)), as.vector(logLik(L_9c$glm_res)), as.vector(logLik(Mod3))),2)
mod_spec = data.frame(rbind(no.obs, aic_mod, log_lik))
rownames(mod_spec) = c("Observations", "Akaike Inf. Crit.", "Log Likelihood")
colnames(mod_spec) = colnames(ttb2)

# Table S1: 
tab2_manuscript = rbind(ttb2, mod_spec)
tab2_manuscript = tab2_manuscript[c("Cash", "Health", "Access" ,"Age", "Male", "High-Educate", "Medium-Educate", "Low-Educate","Employed","Mean-Food","Social Media","District 2", "District 3", "District 4", "District 5","District 6", "Intercept", "Observations", "Akaike Inf. Crit.", "Log Likelihood"),]
stargazer(tab2_manuscript, type = "text", summary = FALSE)
stargazer(tab2_manuscript, type = "latex", summary = FALSE)











